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ABSTRACT.  A number  of  image  analysis  tasks  require  the  registration  of  a 
surface  model  with  an  image.  In  the  case  of  satellite  images,  the  surface 
model  may  be  a map  or  digital  terrain  model  in  the  form  of  surface  elevations 
on  a grid  of  points.  We  develop  here  an  affine  transformation  between  coordin- 
ates of  Multi -Spectral  Scanner  (MSS)  images  produced  by  the  LANDSAT  satellites, 
and  coordinates  of  a system  lying  in  a plane  tangent  to  the  earth's  surface 
near  the  sub-satellite  (Nadir)  point. 
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1.  Introduction. 

Some  image  analysis  tasks  depend  on  the  avail abili 
face  model.  Registration  can  be  accomplished  using  manually  identified  ground 
control  points  or  by  matching  the  real  image  with  a synthetic  image  calculated 
from  the  surface  model  using  assumed  reflectance  properties.  In  either  case, 
the  form  of  the  transformation  from  image  coordinates  to  model  coordinates  must 
be  known.  The  registration  process  is  then  used  to  determine  the  unknown  par- 
ameters of  the  transformation. 

We  show  here  that  in  the  case  of  satellite  images  obtained  by  a mechanical 
scanning  system,  such  as  that  used  on  the  LANDSAT  satellites,  an  affine  trans- 
form applies,  if  small,  second-order  effects  are  neglected.  Such  a transforma- 
tion has  six  parameters  which  depend  on  the  ^♦^ate  of  the  scanning  platform. 

Each  parameter  is  exhibited  as  a function  of  the  components  of  this  state  and 
other  relevant  fixed  quantities.  These  equations  can  then  be  used  to  predict 
transformation  parameters  if  the  state  of  the  scanning  platform  is  known. 

A possible  application  of  automatic  registration  of  images  and  surface 
models  is  the  determination  of  the  parameters  of  a satellite's  orbit.  Unfor- 
tunately, a rigid  body  has  six  degrees  of  freedom  (position  and  attitude)  and 
so  its  state  has  twelve  components  (position,  velocity,  attitude  and  attitude 
rate).  Clearly,  then,  knowing  the  six  parameters  of  the  affine  transformation 
at  one  instant  of  time  is  not  sufficient  to  permit  a calculation  of  the  vehicle's 
state. 

A series  of  determinations  of  the  transformation  for  images  taken  of  dif- 
ferent areas  of  the  earth,  however,  may  permit  determination  of  the  vehicle's 
orbit.  If  we  ignore  small  perturbations,  then  the  position  and  velocity  of 


the  center  of  mass  of  the  vehicle  at  one  Instant  of  time  fully  determine  its 
orbit.  We  prefer  to  use  orbital  parameters  to  describe  this  component  of  the 
state  in  some  circumstances. 
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Basic  Orbital  Geometry^ 


1 


LANDSAT  is  in  a near-polar,  retrograde,  sun-synchronous  orbit  which  is 
nearly  circular.  The  nominal  parameters  of  this  orbit  include  a semi-major 
axis  of  7,294,690  ineters,  that  is,  916,525  meters  above  an  earth  with  equational 
radius  of  6,378,165  meters  and  oblateness  of  1 in  298.3.  The  nominal  period  is 
103.267  minutes,  which  brings  the  sub-satellite  point  back  to  the  same  spot  on 
earth  after  251  orbits  in  18  days.  At  the  equator,  neighboring  sub-satellite 
tracks  are  spaced  159,380  meters.  The  descending  node  is  nominally  passed 
at  9:42  A.M.  The  orbital  inclination  is  nominally  99.092°,  which  brings  the 
satellite  within  e ' 9.092°  of  the  north  pole  at  the  vertex  V (see  Figure  1). 

All  of  the  parameters  drift  with  time  due  to  perturbing  influences  such  as 
solar  wind,  light  pressure,  atmospheric  drag,  non-spherical  distribution  of 
masses  in  the  earth,  effects  of  mass  expulsion,  and  so  on.  The  orbit  is  re- 
adjusted at  time  using  small  gas  discharges  to  maintain  the  positions  of  the 
ground-tracks  within  about  35  km  of  nominal  and  to  prevent  the  time  of  north 
to  south  crossing  of  the  equator  from  drifting  too  far  from  the  nominal  9:42 
A.M.  The  orbital  data  is  derived  from  radio  tracking  information. 

Points  on  the  orbit  may  be  conveniently  referenced  to  the  vertex  V.  The 
orbital  travel  distance  p is  measured  from  it,  and  the  reference  meridian 
passes  through  it.  The  change  in  geographical  longitude  is  measured  from 
the  reference  meridian  (see  Figure  1). 

Ignoring  for  a moment  the  rotation  of  the  earth,  we  find  that  the  nontinal 
position  of  the  satellite,  S^,  lies  at  (geocentric)  latitude  (ji*.  The  nominal 
heading  of  the  satellite  relative  to  the  local  meridian  is  given  by  the  angle 
H^.  The  relationship  between  orbital  parameters  e,  p and  the  geographical  co- 


ordinates  can  be  established  using  products  of  rotation  matrices  [1]. 

Here  we  follow  a more  direct  route  using  spherical  trigonometry. 

Considering  the  right  spherical  triangle  N (see  Figure  1),  applying 

the  sine  theorem,  one  finds  that 


sin  (<j>')/sin  (90° -e)  = sin  (90°-p)/sin  (90°) 


or. 


sin  ' = cos  e cos  p (1 ) 


Next,  from  right  spherical  triangle  V P S^,  one  finds, 
sin  x^/sin  p = sin  H^/sin  e 
and,  applying  the  cosine  theorem  (for  angles). 


cos  X = - cos  cos  90°  + sin  sin  90°  cos  p 
s s s 


or. 


cos  x^  = sin  cos  p 


Hence 
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Equations  (1)  and  (2)  determine  geographical  coordinates  41',  in  tenns  of 
orbital  parameters  e,  p.  Similarly, 

cos  = - cos  cos  90°  + sin  x^  sin  90°  cos  e 


or. 


cos  » sin  x^  cos  c 


Hence, 


tan  = tan  e/sin  p (3) 

Equation  (3)  determines  nominal  heading  in  terms  of  orbital  parameters  c,  p. 

As  it  turns  out,  the  earth  does  rotate  and  so  the  sub-satellite  point  is 
displaced  an  additional  amount  in  the  direction  of  the  local  geographical  para- 
llel, from  point  to  point  S.  The  latitude  remains  unchanged,  of  course, 
while  the  longitude  is  increased  by  x^  and  the  sub-satellite  track  deviates 
by  an  angle  from  the  nominal  direction.  In  order  to  calculate  these  quantities, 
one  must  know  the  angular  velocities  of  the  earth  and  of  the  satellite  in  its 
orbit.  Let  these  quantities  be  and  respectively. 

Since  the  satellite  retraces  its  path  almost  exactly  every  18  days,  after 
completing  251  orbits,  we  know  the  ratio  of  these  two  quantities, 

r * = d X /dp  = 18/251 

w e s e 


(4) 


Then, 


* r p 
e u 


Actually,  is  not  constant,  unless  the  satellite  is  in  a circular  orbit  -- 
we  will  return  to  this  point  later. 

At  latitude  ((>',  the  earth  surface  is  displaced  a distance  cos  (<>'  dt 

in  a time  interval  dt,  where  R is  the  distance  of  the  surface  from  the 
earth's  center. 

The  calculation  of  the  change  in  heading  is  a little  bit  more  complicated. 
If  we  let  the  satellite  heading  H = + H^,  then  one  can  see  that  (Figure  2) 


tan  H = (r  cos  41'  + sin  H^)/cos  H 
lu  S S 


Next, 


tan  Hg  = tan  (H  - H^)  » (tan  H - tan  H^)/(l  + tan  H tan  H^.) 

tan  = r cos  4>'  cos  H,/(l  + r cos  4'  sin  H,)  (7) 

CO)  S (0  S 


Now,  from  the  right  spherical  triangle  P V (Figure  1),  one  obtains 


Sin  (90°  - 4')/sin  (90°0)  = sin  p/sin  H 


That  is. 


cos  ♦'  sin  » sin  e 


4 


Similarly,  one  finds. 
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sin  (90®  - 4i')/sin  (90®)  » sin  p/sin 


That  is, 

cos  = sin  p/sin  (9) 

In  deriving  (3),  we  determined  that  cos  * sin  cos  t,  so  that 

cos  <i'  cos  » sin  p cos  »■  (10) 

Finally,  we  can  use  equations  (8)  and  (10)  to  simplify  the  expression  for 

(7). 

tan  = r cos  t sin  p/(1  + r sin  c)  (11) 

Equation  (11)  determines  in  terms  of  orbital  parameters  p,  p and  the  constant 
r . Note  that  r sin  c is  quite  small  (.0112)  and  can  be  ignored  in  approximate 

U>  (U 

calculations  [1]. 

Also,  now  using  (8)  and  (10),  we  can  simplify  the  expression  for  tan  H (6), 

tan  H " (r  cos‘  + sin  f)/(cos  c sin  p) 


or, 

tan  H ■ [r  (1  - cos^  t cos*  p)  + sin  c]/[cos  r sin  p]  (12) 

U) 
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To  sum  up,  given  e and  <P' ^ we  find  the  orbital  travel  distance  p using 

equation  (1),  the  longitude  relative  to  the  reference  meridian  A = x + x 

s e 

using  equations  (2)  and  (5)  and  the  heading  H * using  equations  (3) 

and  (11 ),  or  (12). 
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Local  Solar  Time  at  the  Point  of  Observation. 

An  immediate  application  of  the  results  developed  so  far  is  the  determina- 
tion of  the  local  solar  time  at  the  sub-satellite  point.  (This  gives  one  some 

idea  of  what  the  position  of  the  sun  is  likely  to  be).  Let  the  time  of  the 

descending  node  be  T^.  That  is,  the  satellite  crosses  the  equator  from  North 
to  South  when  the  local  solar  time  is  T^  (for  LANDSAT  this  is  nominally  9:42  A.M. , 
but  tends  to  vary  as  the  orbit  drifts  and  is  readjusted). 

If  a point  is  observed  when  the  satellite  has  progressed  p in  its  orbit 

from  the  vertex  V,  then  it  still  has  to  travel  through  an  angle  (90°  - p)  be- 

fore reaching  the  equator.  This  will  take  an  amount  of  time  which  can  be  ex- 
pressed in  hours  as  r (90°  - p)/15. 

U) 

Furthermore,  the  point  of  observation  lies  (90°  - \^)  ahead  of  the  point 
of  equator  crossing  in  longitude.  Thus  the  local  solar  time  is  later  by  a time 
which,  when  expressed  in  hours,  comes  to  (90°  - x^)/15.  Finally,  then,  one 
sees  that  the  local  solar  time  at  the  sub-satellite  point  is 

T = T„  + (90°  - X )/15  - r (90°  - p)/15  (13) 

O S U) 

where  tan  x.  = tan  p/sine  (2).  Because  r = 18/251,  one  finds  that  the  first 

S 0) 

term  predominates.  As  a result,  points  North  of  the  equator,  imaged  earl ier 

in  the  orbit,  are  observed  at  local  solar  time  after  T^,  while  points  South 

of  the  equator,  imaged  later  in  the  orbits,  are  observed  at  local  solar  time 

before  T^. 

0 
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4.  The  Scanning  Platform. 

The  satellite  uses  an  oscillating  mirror  to  produce  the  across-the-track 
scan.  Individual  line,  of  the  image  are  obtained  by  this  means.  The  satellite's 
motion  in  orbit  provides  for  the  other  scanning  direction.  Successive  lines 
are  displaced  along  the  sub-satellite  track.  Nominally,  the  optical  system 
points  straight  down  and  the  mirror  scanning  motion  is  perpendicular  to  the 
velocity  vector  of  the  vehicle.  In  practice,  there  are  small  but  significant 
departures  from  this  ideal  state  (Figure  3). 

Pitch  and  roll  are  measured  to  an  accuracy  of  .07"  using  horizon  scanners 
sensitive  to  the  infrared  radiation  (around  14  nm)  emitted  by  the  atmosphere. 

Yaw  is  measur'ed  with  similar  accuracy  using  a gyro  compass.  Pitch  and  roll  are 
maintained  within  .t  .4"  using  the  vehicle's  attitude  contr'ol  system,  while 
yaw  is  maintained  within  * .7°.  A major  component  of  the  attitude  control 
system  is  a set  of  inertia  wheels  which  are  used  in  order  to  keep  down  gas 
expend i ture. 

An  attempt  is  also  made  to  minimize  rates  of  change  of  attitude  which  re- 
sult from  adjustments.  The  maximum  attitude  rates  are  .015  degree/second. 
Attitude  rates  are  estimated  from  time-histories  of  measured  attitudes.  For 
further  information  on  the  scanning  platform  and  its  motion,  see  references 
[1  - 4]. 

Ground  tracking  information  provides  good  ephemeris  data.  However,  since 
a picture  cell  in  the  image  is  only  about  79  meters  by  56  meters,  one  cannot 
expect  the  position  of  the  satellite  to  be  known  accurately  enough  to  predict 
exactly  which  point  of  the  earth  is  imaged.  Similarly,  on-board  horizon 
sensors  permit  a good  determination  to  be  made  of  the  attitude  of  the  satellite 


I 


platform.  Nevertheless,  these  measurements  are  not  accurate  enough  to  permit 
the  direct  calculation  of  the  ground  coordinates  corresponding  to  a particular 
picture  cell.  Errors  of  several  kilometers  may  be  encountered  when  this  is 
attempted  [3]. 

"Precision  processing"  of  satellite  image  information  entails  the  manual 
identification  of  known  ground  control  points  on  each  image  and  the  derivation 
of  a suitable  transformation  based  on  this  information.  So  far,  this  has 
proved  too  expensive  and  LANDSAT  images  are  "bulk  processed",  that  is,  treated 
as  if  the  calculated  position  and  attitude  of  the  satellite  were  exact.  As 
a result,  the  final  photographic  products  may  have  errors  in  translation  of 
several  kilometers.  Fortunately,  non-linear  effects  introduced  by  this  ap- 
proximation are  small. 

One  might  envision  systems  which  automatically  register  image  information 
with  map  or  surface  model  information.  In  such  a system,  one  has  to  model 
the  imaging  operation  so  that  the  registration  process  can  be  used  to  deter- 
mine the  unknown  parameters,  such  as  satel 1 ite  position  and  attitude.  A clear 
understanding  of  the  scanning  process  is  required  to  accomplish  this. 


5.  Fineness  of  the  Scanner  Model. 


A large  variety  of  effects  contribute  to  the  imaging  transformation. 

Amongst  these  are  large  effects  which  must  be  considered,  such  as  the  motion 
of  the  satellite  in  its  orbit  and  the  rotation  of  the  earth  beneath  it.  There 
are  also  smaller  effects  which  have  to  be  judged  individually.  Some  of  these 
produce  non-linear  effects.  Examples  are  panoramic  scan  distortion  (the  mirror 
scans  evenly  in  angle,  not  tangent  of  the  angle),  perspective  projection  (which 
can  be  dealt  with  only  if  a surface  model  is  available)  and  second-order  effects 
of  errors  in  attitude  of  the  spacecraft.  The  relative  importance  of  these  ef- 
fects has  already  been  discussed  by  others  [1  - 4].  The  most  important  criterion 
for  including  an  effect  in  our  model  was  linearity. 

Fortunately,  all  major  components  of  the  image  transformation  turned  out 
to  produce  linear  transformation  of  image  coordinates.  Second  order,  non-linear 
effects  were  neglected,  but  turn  out  to  contribute  errors  which  are  typically 
smaller  than  a picture  cell  In  size.  Compounding  these  linear  transformations 
leads  to  an  overall  affine  transformation  which  is  easy  with  which  to  deal. 

Such  a transformation  has  six  parameters,  which  may  be  found  using  the  registra- 
tion of  the  image  with  some  surface  information  in  the  form  of  a ?nap  or  a digital 
terrain  model. 

The  six  parameters,  as  one  might  expect,  depend  rather  directly  on  the 
position  of  the  satellite  in  its  orbit,  the  attitude  of  the  scanning  platform, 
the  orbital  velocity,  and  the  mirror-sweep  velocity.  It  is  conceivable  that  a 
system  which  automatically  determined  the  parameters  of  the  affine  transformation 
using  a matching  process  of  real  with  synthetic  images  obtained  form  a terrain 
model,  could  also  then  proceed  to  estimate  the  underlying  orbital  data.  A 
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satellite  equipped  with  such  a system  would  be  able  to  determine  its  position 
or  attitude  more  accurately  than  it  might  using  predicted  ephemeris  data  ob- 
tained from  expensive  ground  tracking  efforts. 
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6.  Nominal  Parameters  of  LANDSAT  Imaging  System  (Orbital  Parameters  Drift). 


QrbTt:  Apogee  = 917  km  Perigee  = 898  km 

Inclination  : 99.1°  (Retrograde  orbit) 

Anomalistic  period  : 103.267  minutes 
[That  is,  251  orbits  in  18  days] 

Equatorial  Earth  radius  = 6378  km 

Polar  Earth  radius  = 6357  km 

Equatorial  speed  of  rotation  = 463,8  m/sec 

Average  ground  track  speed  of  satellite  : 6457  m/sec 


Mirror-Scanner:  Frequency  = 13.260  Hertz 

[That  is,  6 lines  are  scanned  every  73.42  msec] 

One  line  every  12.237  msec 

Lines  spaced  by  = 79.0  meters  at  nominal  height 

390  scans  per  image 

That  is,  2340  scan- lines  per  image 

[This  takes  28.63  seconds  and  covers  = 185  km] 

Pixel  Information:  Instantaneous  field  of  view  : 79  m x 79  m 

Mirror  amplitude  = +2.886° 

Total  scan  distance  : 11.545° 

That  is  = 185  km  at  nominal  altitude 

Pixels  per  line  (nominal)  : 3240 

Sampling  interval  = 9.958  psec 

That  is,  about  55.8  - 56.5  m on  the  ground 

Consequently,  : 6.21  radians/second 

Time  to  scan  (six)  lines  (in  parallel)  = 32.238  msec 

Total  Image  Size:  : 2340  x 3240  = 7,581,600  pixels 


I 
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7.  Image  Coordinate  Transformation. 

Let  the  pixels  be  numbered  sequentially  within  each  scan  line  and  let  the 
scan  lines  be  numbered  sequentially.  Then  will  be  the  number  of  a pixel 
counted  from  the  beginning  of  a scan-line,  while  y^  will  be  the  number  of  a line 
counted  from  the  beginning  of  a particular  image.  (Actually,  this  is  arbitrary 
since  the  scanner  does  not  start  or  stop  at  image  boundaries;  the  continuous 
stream  is  segmented  into  images  by  ground  processing).  These  will  be  called 
satellite  coordinates. 

Now  erect  a coordinate  system  in  the  region  of  interest.  First  construct 
a tangent  plane  and  let  the  x-axis  run  in  the  west-to-east  direction,  and  the 
y-axis  in  the  south-to-north  direction.  Now  add  a z-axis  going  vertically  up 
(we  ignore  the  non-spherical  nature  of  the  earth  and  other  such  minor  effects). 
We  will  use  the  notation  for  points  on  the  surface.  The  satellite 

can  also  be  located  in  this  earth  coordinate  system.  At  some  reference  time 
t^,  it  is  at  (xQ»yQ.ZQ)  and  has  attitude  a(roll),  6(pitch),  and  Y(yaw)  (Figure 
3).  The  three  attitude  angles  will  be  assumed  to  be  small. 

At  time  t^,  the  scanner  will  also  be  at  a particular  point  in  its  scan 
of  the  image.  Let  it  be  scanning  the  pixel  in  the  yjQ-th  line  of  the 

image.  If  the  sensor  were  pointing  straight  down  (that  is,  a = 0 and  e = 0), 
it  would  be  imaging  the  sub-satellite  point  (Xp,yp)  (Figure  4). 

At  this  point  we  introduce  a convenient  artifice,  a spherical  earth  fixed 
relative  to  the  orbit  of  the  satellite.  That  is,  a spherical  surface  which 
is  also  sun-synchronous,  rotating  once  a year.  Later  we  will  take  into  ac- 
count the  fact  that  the  earth  rotates  underneath  the  satellite.  We  will  first 
develop  the  coordinate  transformation  for  the  case  of  a fixed  surface  because 
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it  is  easier  to  understand  this  transformation. 

Here  it  is  convenient  to  refer  pixel  locations  to  the  reference  point 
(x^o'^so^  (Figures  4 and  5). 


xi  = (x. 


"so^ 


and 


yi  = (y 


so 


ys) 


(14) 


Let  the  angular  scanning  velocity  produced  by  the  mirror  during  its  linear 
phase  be  0)^^^  (about  6.21  rad/sec)  and  let  t^  (9.958  psec)  be  the  sampling  in- 
terval during  the  scan,  then,  on  a surface  at  a distance  from  the  satellite 
and  perpendicular  to  the  extension  of  the  optical  axis  of  its  scanning  system, 
we  find  a cross-track  scanning  amplitude  X2  as  follows. 


"^2  " ^“m  ^0  ^s^ 

In  the  along-the-track  direction,  the  motion  of  the  satellite  in  its  orbit 
provides  for  the  scanning  and  so. 


yz  = (‘^s  ^ 


(16) 


where  R is  the  distance  of  the  surface  from  the  center  of  the  earth  (about 
6370  km),  while  is  the  angular  velocity  of  the  satellite  in  its  orbit  (about 
1 .014  mill i-rad/sec)  and  t^  is  the  interval  between  successive  scan-lines  (12.237 
mi  1 1 1 -seconds ) . [Actually  six  lines  are  scanned  simultaneously  every  73.42 
milli-seconds.] 

At  this  point,  we  note  that  because  of  possibly  non-zero  yaw,  the  across- 
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track  scanning  may  not  be  perfectly  perpendicular  to  the  along-track  scan. 

This  skewing  effect  can  be  taken  care  of  as  follows  (Figure  6), 

X3  = X2  cos  Y and  ya  “ >2  " *2  sin  y (17) 


We  still  have  to  deal  with  the  effects  of  roll  and  pitch.  For  small  angles, 
these  will  have  the  effect  of  shifting  the  imaged  area  by  an  amount  proportional 
to  the  product  of  the  angles  and  the  distance  to  the  surface  being  imaged. 
Secondary,  non-linear  effects  (such  as  bending  of  the  scanning  line)  will  be 
ignored,  as  will  non-commutativity  of  rotations. 

Thus  the  effects  of  non-zero  roll  and  pitch  can  be  introduced, 

X4  = X3  - a Zjj  and  y4  = ya  - B Zq 

where  z^  is  the  height  of  the  satellite  above  the  surface  as  before.  The  co- 
ordinate system  above  lies  on  the  tangent  place  of  the  (fixed)  sphere.  One 
coordinate  axis  (y)  points  backward  along  the  sub-orbital  track,  while  the 
other  (x)  lies  at  right  angles  to  it.  We  would  prefer  to  work  with  a system 
which  is  aligned  with  local  north.  The  angle  between  the  local  meridian  and 
the  sub-satellite  track  (on  the  fixed  earth)  is  H^.  We  can  rotate  coordinates 
into  a new  system  as  follows  (Figure  7), 

X5  = X14  cos  + y4  sin  (19) 

ys  * -X4  sin  + y4  cos  (20) 


r 
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In  this  new  coordinate  system,  the  y-axis  points  north  and  the  x-axis 
east.  Finally,  we  are  ready  to  introduce  the  rotation  of  the  earth.  It  has 
no  effect  on  the  value  of  y,  of  course,  but  does  introduce  a shift  in  x which 
depends  on  the  time  when  a particular  pixel  is  imaged.  For  a particular  line 
of  the  image,  this  time  can  be  calculated  relative  to  the  time  t^,  when  the 
reference  line  was  imaged.  For  line  number  y^,  this  time  interval  equals 
t^(ys  - y^jj)  * If'  this  time  interval,  the  earth  has  rotated  in  an 

easterly  direction  by  an  amount  which  depends  on  the  latitude. 

^6  “ ^5  + (tg  R cos  t^)  yi  and  yg  = ys  (21) 

where  is  the  angular  rate  of  the  earth  (about  72.722  micro-radians/sec) 
while  (>'  as  before  is  the  (geocentric)  latitude,  and  R the  distance  of  the 
surface  from  the  center  of  the  earth. 

To  obtain  coordinates  in  the  original  system  (x^.y^),  we  must  add  the 
coordinates  of  the  sub-satellite  point  (x^,yQ), 


X 


e 


^6  + Xq 


ye 


and 


(22) 


1 
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8.  The  Overall  Transformation. 

All  the  partial  transformations  can  now  be  combined, 

Xg  = Xi,  cos  sin  R t^)  cos  yj  + x^ 

Yg  = -X4  sin  + y4  cos  + y^ 

Or. 

Xg  = X3  cos  + ys  sin  + (wg  R t^)  cos  y^  - 
(a  cos  H + 0 sin  H ) z + x 

Yg  = sin  + yj  cos  - (-a  sin  + 0 cos  H^)  z^  + y^ 

That  is, 

Xg  = cos  (H^  + y)  X2  + sin  yg  + (wg  R t^^)  cos  4''  yi  + - 

(a  cos  + 0 sin  Hg)  Zg 

Yg  = -sin  (Hj  + y)  X2  + cos  Y2  + Yq  - (-a  sin  + 0 cos  H^)  z^ 

Or, 

Xg  = cos  (Hj  + y)U^  Zq  tg)xi  + [sin  R t^)  + cos  ♦'(ugR  t^^)]  yj  + 
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9.  Form  of  the  Transformation. 


The  transformation  is  of  the  form. 


= a Xi  + b yi  + c 

(23) 

= d Xj  + e yi  + f 

(24) 

This  is  an  affine  transformation,  where  the  six  parameters  are  given  by 


a = 2^  tg)  cos  (H^  + y)  (25) 
b = (ws  R t^)  sin  + (wg  R t^^)  cos  (26) 
c = Xq  - (a  cos  + 6 sin  H^)  (27) 
d = -(co^  2^  t^)  sin  (H^  + y)  (28) 
e = (ws  R t^)  cos  (29) 
^ = Yq  -(-“  sin  + B cos  H^)  (30) 


We  can  use  these  equations  to  predict  approximate  transformation  parameters  from 
estimated  values  of  satellite  position,  attitude  and  velocity  in  orbit.  Con- 
versely, if  we  can  use  ground  control  points  or  digital  terrain  models  to 
determine  the  coefficients  of  the  transformation  more  precisely,  we  can  try  and 
improve  the  estimates  we  have  of  satellite  position  and  attitude. 

From  the  form  of  the  equations  for  c and  f it  becomes  imnediately  clear, 
however,  that  there  are  some  limits  to  this  process.  That  is,  one  cannot  dis- 
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tinguish  in  our  model  between  displacements  of  the  satellite  across  the  track 
and  small  roll  errors.  Similarly,  displacements  along  the  track  have  the 
same  effects  as  small  pitch  errors.  Thus  two  of  the  six  components  of  position 
and  attitude  cannot  be  found  this  way. 
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10.  Attitude  Rates. 

Pitch,  roll  and  yaw  drift  during  the  scanning  of  a single  image.  The 
rates  are  less  than  .015  degrees/second.  A constant  rate  of  change,  s of  pitch 
has  the  same  effect  as  a change  in  along-track  ground  velocity  of  6 z^.  That 
is,  equation  (16)  becomes, 

ya  = (o)^  R + e z^)  t^  yi  (31) 

The  transformation  is  altered  only  in  the  appearance  of  (u^  R + B z^)  t^  in 
place  of  (u)^  R t^).  Typical  values  for  w^R  and  Bz^  are  6458  and  32  meters/ 
second  respectively  (when  B = .002  degrees/second).  This,  then,  is  a small 
but  noticeable  change  {'  .5%). 

A constant  rate  of  change,  a,  of  roll  has  little  effect  on  the  scanning 

of  a single  line  since  w z >>  a z . Successive  lines,  however,  are  shifted 

mo  0 

laterally  by  a z^  t^.  That  is,  equation  (17)  becomes, 

X3  = X2  cos  Y + (a  Z^  t^)yi  and  ya  = yz  - X2  sin  y (32) 

This  causes  additional  skewing  of  the  image. 

Here,  typical  values  are  u z t - 56.5  meters  and  a t^  t.  ' .4  meter 

m 0 s ox. 

(for  a = .002  degree/second).  Again  we  see  a small,  but  noticeable  change 
('  .7%),  resulting  in  additional  skewing  (-  .4%). 
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Finally,  the  transform  parameters  are  now: 


a = (u)^  2^  t^)  cos  (H^  + y)  (33) 

b = (u)  R + 0 z ) t sin  H + (a  z t ) cos  H + (u  R t ) cos  <(1 ' (34) 

S 0^  S Oik  S6iv 

c = X Q -(a  cos  + e sin  H^)  z^  (35) 

d = -(tOjj^  t^)  sin  (H^  + y)  (36) 

e = (<D^  R + g Zg)  tj^  cos  - (o  z^  tj^)  sin  (37) 

= Yq  -(-a  sin  + B cos  H^)  z^  (38) 


We  see  here  that  the  transformation  parameters  depend  on  the  timing  (t^.t^) 
of  the  scanning  system,  the  mirror  scan  velocity  (t^^).  the  position  of  the 
satellite  relative  to  the  tangent  plane  coordinate  system  (^Q’yQ’^p),  the 
orbital  velocity  vector  (as  it  affects  and  H^),  the  attitude  of  the  scanner 
(a,  g,  y).  the  attitude  rates  (a,  g,  y).  and  the  latitude,  41'. 

The  vertical  component  of  the  velocity  vector  (altitude  rate)  and  also 
the  yaw  rate  do  not  contribute  to  linear  changes  in  the  imaging  transformation. 
The  first  produces  a change  of  lateral  scale  from  one  end  of  the  image  to  the 
other,  the  second  a tilt  of  image  lines  at  one  end  of  the  image  relative  to 
those  at  the  other  end.  Such  small,  non-linear  effects  are  ignored.  Except 
for  these  two,  however,  all  twelve  components  of  the  state  of  the  scanner 
platform  influence  the  transformation  parameters. 

If  yaw  rate  and  altitude  rate  jre  included,  one  finds  small  terms  in  xjy, 
(and  yp.  The  transformation  is  then  no  longer  an  affine  transformation. 


For  small  regions,  the  effects  of  these  terms  can  be  ignored  --  for  images 
which  are  large  fractions  of  a standard  LANDSAT  image,  they  can  not.  In  the 
latter  case,  one  has  to  include  other  non-linear  terms  we  have  ignored  in  any 
case  and  then  the  transformation  can  be  expressed  with  sufficient  accuracy  by 
two  second-order  polynomials. 
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11.  Figure  of  the  Earth. 

I' 

To  calculate  the  displacement  of  points  due  to  the  rotation  of  the  earth, 
and  to  relate  the  geocentric  distance  of  the  satellite  to  its  altitude  above 
the  surface,  one  needs  to  be  able  to  calculate  the  distance  of  a point  on  the 
earth's  surface  from  the  earth's  center.  To  a first  approximation,  a meridonal 
cross-section  through  the  earth  is  an  ellipse  (Figure  8)  with  semi-major  axis, 
a = 6,378,165  meters  at  the  equator;  and  semi-minor  axis,  b = 6,356,783  meters 
at  the  poles. 

If  we  introduce  a coordinate  system  with  the  x-axis  along  the  semi -major 
axis  and  the  y-axis  along  the  semi-minor  axis,  then  the  geocentric  latitude, 

<j('  is  defined  by, 

tan  ({>'  = y/x  (39) 

The  more  commonly  used  geographic  latitude,  <t>,  is  the  angle  between  a local 
normal  to  the  surface  and  the  equatorial  plane.  Thus, 

-1/tan  4>  = ^ (40) 

i Using  the  equation  for  the  ellipse, 

1 

I {x/a)2  + (y/b)2  = 1 


one  finds  that 


T * 


r 

I 
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SO  that 


tan  ,>*  = (b/a)2  tan  4> 


(41) 


Further, 


X = (- 


ab 


tan^  41' 


y = ( 


ab 


+ a^  tan'-^  .ji' 


-)  tan  1(1' 


The  distance  of  a point  from  the  center,  then,  is. 


R = /x-  +~1P  = 


ab 


1^2] 


The  height  of  a surface  feature  above  mean  sea-level  must  be  added  to  this 
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12»  Orbital  Velocity. 

Since  one  scanning  motion  depends  on  the  satellite's  angular  velocity  in 
its  orbit,  it  is  useful  to  relate  this  to  other  quantities.  Using  Kepler's 
second  law,  one  immediately  sees  that 


r2 


dt 


h 


(43) 


where  r is  the  radius  vector  (geocentric  altitude  of  the  satellite),  while  h 
is  a constant.  Now  the  area  of  an  ellipse  with  semi -major  axis,  a,  and  semi- 
minor  axis,  b,  is  irab,  so  that 


where  T is  the  complete  period  of  the  satellite.  Using  for  the  angular 
velocity  of  the  satellite,  one  finds. 


_ 2Tr  ab  _ . ab 
- T?7  ~ 


(44) 


where  the  average  angular  rate  h = 2v/J.  Further, 

h2a3  = li  = GM  (45) 

by  Kepler's  third  law,  where  G is  the  gravitational  constant  and  M is  the  mass 
of  the  earth,  p : 396.08  x 10’2  radian2  meter3/second2. 


I 


k 


E ■ 
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For  orbits  with  small  eccentricity,  there  is  very  little  difference  be- 
tween a and  b,  so  perigee,  r , and  apogee,  r^,  are  given  instead,  where 

P ® 


rp  = a (1  - e) 
r^  = a (1  + e) 


and 


a2  (1  - q2)  = b2 


Thus, 


a = (r,  * rp)/2 


a p 


(46) 


(47) 


(48) 


Note  that  perigee  and  apogee  are  frequently  given  as  distances  from  the  surface 

of  the  earth  and  the  equatorial  radius  of  the  earth  (6,378,165  m)  has  to  be 

added  to  this  in  order  to  find  r and  r . For  a more  detailed  analysis,  see 

P « 

Lyndanne’s  modification  of  Brouwer's  analysis  of  satellite  orbits. 


LL 
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13.  Maj)  Distortion. 

So  far  we  have  ignored  the  fact  that  the  spherical  surface  of  the  earth 
cannot  be  represented  on  a planar  map  without  distortion.  Up  to  this  point, 
coordinates  have  been  referred  to  a hypothetical  plane  tangent  to  the  earth's 
surface  at  a point  in  the  region  of  interest.  Typically,  a digital  terrain 
model  will  be  derived  from  a map  with  a different  projection  and  a transformation 
must  be  established  between  the  two  coordinate  systems. 

It  can  be  shown  that  for  typical  map  projections  such  as  transverse 
Mercator  and  conformal  oblique  axis  cylindrical  there  exists  a small  rotation 
of  map  coordinates,  where 

sin  H : sin  <ti  sin  (e  - e„)  (49) 

moo 


Here  the  projection  is  centered  on  a point  at  longitude  and  latitude  41^ 
and  the  point  of  interest  is  at  longitude  e and  latitude  41.  Consequently, 
the  map  coordinates,  (Xjj^,yjj^),  are  related  to  the  geographical  coordinates  on 
the  tangent  plane  (Xg,yg)  by 


cos 

«m 

-sin  H_ 
m 

. 1 

s 

sin 

cos 

•'m 

m 

m 

(50) 


where  s the  scale  of  the  given  map.  There  will  also  be  a small  scale  change 
which  varies  with  1/cos  (0  - 0^)  for  transverse  Mercator  and  with  1/cos  ((j*  - 4"^) 
for  conformal  oblique  axis  cylindrical  projection.  Typically,  this  effect  is 


V 
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so  small  that  it  may  be  Ignored. 

The  affine  transformation,  (23)  and  (24)  must  be  pre -multiplied  by  an 
augmented  rotation  matrix  M to  correctly  relate  satellite  image  coordinates 
to  map  coordinates. 


cos  -sin  0 

m m 

sin  cos  0 

m m 

0 0 1 


(51) 
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14.  Example  of  Map  Rotation. 

The  national  map  of  Switzerland  is  based  on  a conformal  oblique  axis 
cylindrical  projection  with  center  at  the  city  of  Bern. 

0 : 7‘’26'20"  A : 46“57'10" 

0 ^0 

The  region  of  interest  covered  by  an  available  terrain  model  lies  at 

0 : 7°8'  ^ : 46‘'15' 

The  map  rotation  then  is  -.225°  (-.00393  radians).  The  transformation  matrix 
becomes 

1.0  +.00393 

-.00393  1.0 

(The  scale  error  is  less  than  one  part  in  10,000  and  can  be  ignored).  For  more 
details  regarding  suitable  map  transformations  see  Col vocoresses  paper  on  the 
"Space  Oblique  Mercator"  projection. 


! 


-34- 


I 


15.  Numerical  Example  --  Estiiratlnq  Transform  Parameters. 

LANDSAT  image  number  1078-09555,  produced  1972/0ctober/9  shows  a region 
of  Switzerland  including  a mountain  range  called  "Dent  de  Hordes".  The  image 
annotation  data  suggest  that  at  the  nadir  the  geographic  latitude  was  45.9197“ 
and  the  heading  193.11324°.  Using  equation  (41),  we  see  that  the  geocentric 
latitude  of  the  nadir  point  is  45.7274°  and  so,  using  equation  (8),  it  appears 
that  the  orbital  inclination  must  have  been  c = 9.11°. 

The  altitude  is  given  as  915,724  meters  near  the  region  of  interest,  which 
lies  at  an  average  of  1700  meters  above  average  sea  level.  So  z^  : 914  km.  The 
angular  velocity  of  the  satellite  is  slightly  above  its  average  rate, 

“s  “ 24  X 60  X ^ ^ ^ ^ 1.00967  radians/second 

The  region  of  interest  lies  at  geographic  latitude  (|)  = 46.25°  (and  thus  at 
geocentric  latitude  41'  = 46.06°),  while  the  scanner  at  that  time  is  above 
geographic  latitude  41  = 46.40°  (that  is,  geocentric  latitude  41'  = 46.21°). 

Further, 

“e  = 24  X lo  X 60  i^^dians/second 
u)„  = 6.21  radians/second 


tg  = 91958  psecond 


I 
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The  image  annotation  also  gives. 


a (roll)  = -.20370“ 

6 (pitch)  = .06688“ 

Y (yaw)  = .23387“ 

a = -.00160“/second 

B = -.00109“/second 

Y =.00189“/second 

The  earth  radius  (equation  42)  near  the  region  of  interest  is  6,367,081  meters. 

Adding  the  average  elevation  above  sea  level,  we  get 

R = 6,368,800  meters 

Using  equation  (1),  one  finds  that  p = 43.021“  and  by  equation  (8),  that  = I 

13.226“.  Further,  (H^  + y)  = 13.460“.  Then  also,  0)^^^  t^  = 56.521  meters, 

(u)  R + 6 Zn)f|,  = 79.582  meters,  a z t = -.312  meter,  and  w R t = 5.667 
meters. 

So,  finally, 

a = 54.969  b = 21.837  c = x + 2919 

0 

d = -13.156  e = 77.543  f = y - 1782 

0 

See  Figure  9 for  a graphical  illustration  of  this  transformation.  It  shows  as 
a parallelogram  the  region  of  the  surface  scanned  when  an  image  with  a number 


of  lines  equal  to  the  number  of  pixels  per  line  is  gathered. 

We  can  take  the  inverse  of  the  2x2  sub-matrix  and  obtain. 
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Xi  = .01704  - .00480  + x' 

* e e 0 

y,  = .00289  x„  + .01208  y„  + y' 

* € 6 0 

where. 


See  Figure  10  for  a graphical  illustration  of  this  transformation.  It 

shows  as  a parallelogram  the  appearance  in  the  image  of  a square  region  on  the 

( 

ground  aligned  v-^th  the  north-south  and  east-west  axes.  Finally,  if  we  have 
a terrain  model  on  a 100  x 100  meter  grid,  such  that  x^  = 100  * i and  y^  = 

100  * j,  then, 

Xi  = 1.704  i - .480  j + x' 

0 

yi  = .289  i - 1.208  j + y^ 

■r 

Finally,  we  have  to  introduce  the  map  distortion  by  post-multiplying  by  the 
inverse  of  the  map  transformation  matrix  introduced  earlier  (the  inverse  of  a 
rotation  matrix  equals  its  transpose).  The  matrix  then  becomes 


1.702 

00 

1 

.294 

1.207 
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16.  Detennininq  Orbit  Parameters  from  Transfonnatlon. 

Many  parameters  appear  in  the  equations  for  the  six  coefficients  of  the 
affine  transformation.  Some  are  known  accurately,  others  only  approximately. 
For  example,  t^  and  t^^  are  fixed  fairly  accurately  by  electronic  oscillations 
on  board  the  satellite.  The  angular  rate  of  the  earth  is  constant  and  ac- 
curately known,  while  the  angular  rate  of  the  satellite,  depends  on  its 
altitude,  semi-major  axis,  orbit  eccentricity,  and  orbital  period.  The  angular 
rate  of  the  scanning  mirror,  varies  somewhat  during  each  scan,  though  in 
a fairly  repeatable  fashion.  The  earth  radius,  R,  can  be  calculated  with 
sufficient  accuracy  if  the  latitude,  <J)',  is  known  approximately. 

It  is  reasonable,  then,  to  assume  that  one  can  accurately  predict  a value 

for  a)„  R t„  cos  4i'.  Then  let 
e t 

b'  = b - 0)  R t„  cos  <j)' 
e t 

The  following  equations  permit  the  determination  of  useful  satellite  parameters 
from  transformation  parameters: 

+ Y = tan'i(-cl/a)  (52) 

((i)_  z t ) = /a^  + d^  (53) 

m 0 s ' 

/(o)^  R + 6 Zq)^  + (o  = /(b’  )^  + e^  (54) 


If  the  heading  can  be  estimated,  one  finds 
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(Wj  R + 6 Zq)  = b'  sin  + e cos 


(&  Zq)  = b'  cos  sin 


(55) 

(56) 


Alternatively,  the  heading  can  be  calculated  if  the  attitude  rates  are  known, 


_ b'  (03  R + B z)  - e (aZ) 
“s  = B-'-Ta  zj  4 e{u>3  R > ft  z^) 


(57) 


Clearly  there  are  limitations  to  this,  since  the  full  state  of  the  scanning 
platform  cannot  be  ascertained  from  six  parameters  alone.  A series  of  such 
measurements  is  needed  to  determine  all  twelve  components  of  state. 
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i 
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17. 


Numerical  Example. 


Suppose  the  transformation  matrix  applicable  to  a 100  x 100  meter  grid 
was  found  by  image  registration  techniques  using  synthetic  images, 


1.694 

-.512 

.300 

1.216 

for  any  area  with  a map  distortion  as  defined  previously.  Post-multiplying 
by  the  map  transformation  matrix  gives 

1.693  -.505 

.295  1.217 

The  area  is  at  geographical  latitude  (p  = 46.25°.  The  geocentric  latitude  is 
then  ((i*  = 46.06°  and  so  R = 6,367,081  meters.  If  the  region  of  interest  lies 
at  an  average  altitude  of  1700  meters,  and  t^^  = 1/(6  x 13.62)  seconds,  then 

“e  ^ ~ meters/scan-line 

The  inverse  transformation  matrix  then  is 

55.08  22.86,  \ 

I 


-13.35 


76.63 
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From  this  one  finds, 

Hj.  + Y - 13.624“ 

“m  ^0  ^s  " 56.67  meters/pixel 

And,  if  the  attitude  rates  are  assumed  to  be  very  small, 

- 13.874“ 

and  so. 


(“•5  Rtj)  ■ 78.93  meters/scan-line. 


_41. 

18.  Using  Ground  Control  Points. 

1 

An  alternate  method  of  estimating  transformation  parameters  is  based  on  i 

'1 

the  identification  of  points  of  known  ground  position  in  the  image.  Since  | 

I 

the  affine  transformation  has  six  parameters,  one  needs  to  locate  three  such 
points.  Let  the  coordinates  of  the  points  be  (xj,  yj),  (x2,  yj)  and  (X3,  y3) 
in  the  image  and  (xj,  yj),  (x2,  ya)  and  (X3,  y^)  on  the  map,  then 


x!  = a x^  + b y^  + c 
y*  = d x^  + e y^  + f 
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If  more  than  three  points  can  be  identified,  better  accuracy  is  available 
using  a least-squares  procedure.  That  is,  if 


xi  yi  1 

X2  yz  1 


Then  a good  set  of  values  for  the  parameters  is 


a 

x{ 

b 

c 

d 

y{ 

e 

yz 

f 

y’ 

•^n 

Typically,  the  accuracy  obtained  by  fitting  a discrete  set  of  ground  control 
points  has  been  found  to  be  inferior  to  the  area-based  matching  of  real  and 
synthetic  images. 
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